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We describe the use of the Density Matrix Renormalization Group method as a means of ap- 
proximately solving large-scale nuclear shell-model problems. We focus on an angular-momentum- 
conserving variant of the method and report test results for the nucleus 48 Cr. The calculation is able 
to reproduce both the ground state energy and the energy of the first excited state, by diagonalizing 
matrices much smaller than those of the full shell model. 
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I. INTRODUCTION 

In the traditional nuclear shell model [J, the low- 
energy structure of a given nucleus is described by as- 
suming an inert doubly-magic core and then diagonaliz- 
ing the effective nuclear hamiltonian within an active va- 
lence space consisting of at most a few major shells. De- 
spite the enormous truncation inherent in this approach, 
the shell-model method as just described can still only be 
applied in very limited nuclear regimes, namely for those 
nuclei with a sufficiently small number of active nuclc- 
ons or a relatively low degeneracy of the valence shells 
that are retained. For heavier nuclei or nuclei farther 
from closed shells, one must truncate further to reduce 
the number of shell-model configurations to a manage- 
able size. 

An attractive truncation possibility is provided by 
the Density Matrix Renormalization Group (DMRG), a 
method initially developed for low-dimensional quantum 
lattices Q. Subsequently, the method was extended to 
finite Fermi systems and applied to ultrasmall supercon- 
ducting grains to problems in quantum chemistry 
and to two-dimensional electrons in strong magnetic 
fields p| . The great success of these applications suggests 
that it could also prove useful for obtaining accurate ap- 
proximate solutions to the nuclear shell model, in cases 
where exact diagonalization is not feasible. 

The DMRG method involves a systematic inclusion of 
the degrees of freedom of the problem. When treating 
quantum lattices, real-space sites are added iteratively. 
In finite Fermi systems, these sites are replaced by single- 
particle levels. At each stage, the system [referred to as 
a block] is enlarged to include the additional site or level. 
This enlarged block is then coupled to the rest of the 
system (the medium) giving rise to the superblock. For 
a given eigenstate of the superblock (often the ground 
state) or perhaps for a group of important eigenstates, 
the reduced density matrix of the enlarged block in the 
presence of the medium is constructed and diagonalized 
and those states with the largest eigenvalues are retained. 
This method of truncation is guaranteed to be optimal 



in the sense that it maximizes the overlap of the ap- 
proximate (truncated) wave function with the targeted 
superblock wave function prior to truncation. 

The earliest application of the DMRG method in the 
context of nuclear structure was based on the use of a 
particle- hole variant of the method [f| . In this approach, 
the group of single-particle levels was divided into four 
subgroups, neutron particle, neutron hole, proton par- 
ticle and proton hole. The results of those calculations 
were not terribly encouraging, however, as they required 
that a substantial fraction of the full shell model space 
be retained for a reasonable reproduction of the results 
from exact diagonalization. 

More recently a more traditional DMRG study was 
reported that did not adopt the particle-hole strategy • 
The method worked extremely well for the nucleus 28 Si, 
giving an accurate reproduction of the exact shell-model 
results with a significantly reduced basis. In contrast, the 
results for 5e Ni were not nearly as good, as the method 
seemed to be converging to a solution energetically still 
far from the exact ground state. 

One of the limitations of the above studies is that they 
did not preserve symmetries throughout the iterative in- 
clusion of single-particle levels. Since the density matrix 
procedure involves a truncation at each of the iterative 
stages, there is a potential to lose these symmetries and 
the associated correlations. In such cases, the iterative 
procedure must incorporate not only dynamical corre- 
lations but also the kinematical correlations associated 
with symmetry restoration. 

This point has been recognized for some time in DMRG 
applications outside the nuclear domain. McCulloch and 
Gulacsi [8j discussed how to reformulate the DMRG to 
maintain symmetries throughout. They referred to this 
as the non-Abelian DMRG method and showed that it 
could produce more accurate results than the traditional 
Abelian method while using smaller spaces. 

In this work, we adopt a strategy whereby angular 
momentum is preserved throughout the iterative DMRG 
process. We refer to this as the JDMRG method. It was 
briefly described as part of a recent review Q. Subse- 
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quently, the method was applied in nuclear physics for 
the first time in the context of the Gamow Shell Model 

An outline of the paper is as follows. In Sec. II, we pro- 
vide a brief overview of the traditional Abelian DMRG 
method and then in Sec. Ill describe the modifications 
required to maintain rotational invariance as an exact 
symmetry. In Sec. IV, we briefly review the test problem 
we treat and describe in some detail our implementation 
of the JDMRG method for this problem. Section V de- 
scribes our results. Section VI summarizes the principal 
conclusions reached in this study and outlines directions 
for future work. 



II. OVERVIEW OF THE DMRG METHOD 

We begin by describing the traditional Abelian DMRG 
method, as has been used in almost all applications to 
date. We focus on finite Fermi systems, for which the 
degrees of freedom being systematically incorporated by 
the method are single-particle levels. 

Assume we have treated a set of L levels and that the 
total number of states we have kept to describe that block 
is m. 

We now add the next level, the L+ 1 st , which admits s 
states. In the traditional DMRG approach, the states we 
build in the enlarged block are simply products of those 
in the original block and in the extra level. Thus, when 
we enlarge the block to include the L + 1 st level, it results 
in a new block with m x s states. 

Renormalization Group methods, whether Wilson's 
original numerical algorithm or White's Density Matrix 
version, then implement a truncation of the states of the 
enlarged block to the same number m as before the block 
enlargement. 

In Wilson's numerical RG procedure, the truncation 
involves diagonalizing the hamiltonian in the m x s- 
dimensional space of the enlarged block and retaining 
its lowest m eigenstates. 

In White's DMRG approach, the truncation is imple- 
mented through a very different procedure, schematically 
illustrated in figure 1. Consider the enlarged block B' 
(the block B plus the added level) in the presence of a 
medium M that approximates the rest of the system. 
White's truncation is carried out based on the impor- 
tance of the block states in a selected set of target states 
of the full superblock, i.e., the eigenstates of B' coupled 
to M. 

Assume here that we target only the ground state of 
the superblock in the truncation, 

!*»«)= E E %K>BbV, (i) 

mxs j — 1, t 

obtained by constructing and then diagonalizing the su- 
perblock hamiltonian. Here refers to a state in the 
block and \j)m to a state in the medium. Furthermore, t 
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FIG. 1: Schematic illustration of the DMRG growth proce- 
dure, in which a block B is enlarged to B' in the presence of 
a medium M. 

denotes the number of states of the medium. If we con- 
struct the ground-state density matrix for the enlarged 
block, 

($> = E > ( 2 ) 

3=1, * 

diagonalize it, and truncate to the m states with the 
largest density- matrix eigenvalues, we are guaranteed 
to achieve the optimal approximation to the superblock 
ground state. 

To target a group of states, we would construct a mixed 
density matrix containing information on the block con- 
tent of all of them. For example, if we wished to target 
both the ground state and the first excited state, weight- 
ing them equally, we would construct 

& = \ ( E *«*<'i + E <|, v" '. I > ( 3 ) 

where ^ as above refers to the ground state wave function 
and $ correspondingly to the wave function of the first 
excited state. 

Truncation to the most important m states in the en- 
larged block is accompanied by the renormalization of 
all operators for use in the truncated space. If the eigen- 
states | a) associated with the m lowest density matrix 
eigenvalues are 

\a) = c i\^ > a = l,...,m, (4) 

i—l, mxs 

then the matrix elements of any operator O must be 
renormalized according to 

<a|0|/3)=5>rcf <«>. (5) 

A key step in the growth procedure is to calculate 
at each step the matrix elements of all necessary sub- 
operators of the hamiltonian and to store them. This 
includes all one- and two-point operators, 

a\, a\aj, atoj, + h.c, (6) 
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FIG. 2: Schematic illustration of the ordering of levels (or 
sites) in a chain to be systematically included in the DMRG 
growth procedure. 

and all three- and four-point operators required to build 
the hamiltonian matrix, 

= 4 v Hkl a \ a ] a l + h - c - > 

ijl 



O 1 = | ^2 Vi 3 kl 4 a ] a l a k , (?) 
ijkl 

where Vijki — (ij\V\kl). Having this information for the 
block and the additional level enables it to be calculated 
for the enlarged block. 

As an example, consider the one-body operator a^ap. 
Its matrix elements in the enlarged product space \i,j) — 
\i)\j), where |z) belongs to block B\ and \j) belongs to 
5 2 , are given by 

(hj\a f a at3\k,l) = (i\a) a ap\k) S jt 

+ (j\ a l a l3\ l ) $ik 

+ (-)^<ijo+|fc> V\ap\l) 

{-Y*(i\a p \k) (j\al\l) , (8) 

where n, is the number of particles in state Note 
that the matrix elements of a two-point operator in the 
product space is expressed in terms of matrix elements 
of one- and two-point operators in the spaces of the two 
component blocks. 

We now have the basic tools in hand to discuss how the 
traditional DMRG procedure is implemented. It involves 
a series of steps. The first is to choose an active space 
and a hamiltonian to act in that space. Next we choose 
an order in which the single-particle levels are included. 
Typically, this is done by organizing the levels in a chain, 
as schematically illustrated in hgurc 2. 

The following step is usually referred to as the warmup 
phase. This step involves making an initial guess for the 
optimal structure in each size block for a given choice 
of the number of retained states to. For the eight levels 
shown in figure 2, we would make a first guess at the 
structure of the blocks 5(1,2), 5(1,3), 5(1,4), 5(1,5), 



5(1,6), 5(1,7) and 5(1,8). Having a guess for the 
optimal to states in each block, we can calculate (and 
store) all the necessary matrix elements in those trun- 
cated blocks. The warmup guess can be implemented in 
a variety of ways, depending on the specific problem. 

The next step, and indeed the heart of the DMRG 
method, is the sweep phase, which we now describe by 
focusing again on the chain of eight levels in figure 2. 

We begin the sweep phase by growing level 8 into a two- 
level block 5(7,8). To determine the optimal m states 
in that enlarged block, we immerse it in a medium con- 
sisting of the block 5(1,6) from the warmup phase. We 
diagonalize the superblock hamiltonian for the full sys- 
tem, then calculate the reduced density matrix associated 
with the enlarged block, and renormalize all operators to 
act in the truncated space associated with the to largest 
eigenvalues. 

We then continue the process, enlarging to 5(6, 8), im- 
mersing it in a medium consisting of 5(1,5) and then 
truncating and renormalizing 5(6,8) based on its re- 
duced density matrix eigenvalues. 

This process is continued until all size blocks have been 
treated. 

At this point, we simply reverse the direction of the 
sweep, now growing the system from left to right, always 
treating the enlarged block in a medium consisting of 
the rest of the levels as obtained in the previous sweep. 
This process is continued until acceptably small changes 
from one sweep to the next are achieved. At this point, 
we have hopefully arrived at an optimal block structure 
and an optimal description of the system for the assumed 
value of to. 

The calculation is then carried out as a function of to, 
until acceptably small changes with increasing values are 
achieved. At this point, we hopefully have a good approx- 
imation to the exact solution of the problem. Obviously, 
the solution is especially useful if this good approxima- 
tion can be achieved while retaining only a small fraction 
of the full Hilbert space. 



III. THE JDMRG METHOD 

We now describe the modifications to the traditional 
DMRG algorithm required to preserve rotational symme- 
try throughout. 

As noted earlier, the DMRG systematically grows the 
system by adding additional degrees of freedom. This 
is most conveniently accomplished by adding doubly- 
degenerate levels, either Nilsson-like levels or the +to 
and —to levels of an nlj multiplet. A great advantage of 
this is that every level looks very similar to every other, 
greatly simplifying the algorithm by which they are com- 
putationally included. 

In the JDMRG, it is essential that throughout the 
procedure the states have definite angular momentum. 
To accomplish this, it is critical that we add at each 
stage states of definite angular momentum. This is most 
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readily accomplished by always adding a complete shell, 
rather than a partial shell or a non-spherical single- 
particle level. 

The goal of a calculation is to construct and diago- 
nalize the hamiltonian in the the space in which all lev- 
els have been included. As noted earlier, the traditional 
DMRG works in a product (or m-scheme) representation. 
In this case, we need to build the matrix elements of all 
hamiltonian sub-operators directly in the m-scheme, as 
this is what is needed to build the matrix elements in 
the enlarged system and ultimately the matrix elements 
of the superblock hamiltonian. For each size block those 
m-scheme matrix elements must then be stored, for use 
in the following sweep. 

In an angular-momcntum-conserving approach we 
must instead calculate and store throughout the itera- 
tive process the reduced matrix elements of all relevant 
sub-operators of the hamiltonian, namely 



(a'K', 0'L'(J')\\ [a%] x \ \aK , (3L(J)) = 

c c ( ^K'+L' + J+X <j<fi / K' J 1 L' 

Opp'dLU (-) \ J K X 

x(a'K'\\ (a\a j y\\aK) + 

■^K'+L + J'+X J-ji 



+ 



x (f?L'\\ (ofo) \\/3L) 

( K' V J 1 
+ (-)"-* JJ'\\ K L J 
[ * J A 

x (a'K'\\ a\ \\aK) x {P'L'\\a 3 \\0L) + 

( K' V J' 

+ (-y- 3+X {-) n " K JJ'Xl K L J 

{ 3 i A 

x (a'K'WajWaK) x (/?'£' || a\ \\0L) , 



V J' K' 
J L A 



(10) 



aj, [a\~ aj ] K , [a]a]] K , 



6} = 



ijkK 



(1 + 6^(1 + S k i) KV« kl 
(Ha]] K ~a k ) 1 , 



ijklK 



x ([a\a)] K [a fe a ; ]^)° , (9) 



and the relevant hermitean conjugates. Here, K 
\/1K + 1, a jm = (-) j ~ m aj- m and 



where n a K is the number of particles in state \aK) . 
Knowing the reduced matrix elements of the one- and 



two-point operators a\ , [a J a j ] , [a] a] ] in the two com- 
ponent blocks enables us to determine the reduced matrix 
element of interest for the enlarged block. 

Otherwise, the JDMRG involves the same set of basic 
steps as the traditional DMRG. 



IV. THE TEST MODEL 

The test results we report in this paper are for the 
nucleus 48 Cr, treated as four valence neutrons and four 
valence protons outside doubly-magic 40 Ca. We restrict 
both the valence neutrons and the valence protons to the 
orbits of the l/2p shell. Furthermore, we use in these 
calculations the shell model hamiltonian KB3, for which 
exact results are available. The size of the full space for 
48 CV is 1,963,461 states. Of these, 41,355 are 0+ states, 
182,421 are 2+ states, etc. 



V%h = (ij(K)\V\kl(K)). 



Some specifics of our implementation for this 
test problem 



As an illustration of how this proceeds, consider the re- 
duced matrix elements of the coupled one-body operator 
[ajaj]^ in an enlarged space obtained by coupling blocks 
B\ and B 2 . Denoting the coupled states as \aK, (3L( J)), 
where \aK) refers to the state in B\, \0L) to the state in 
B2, and J to the coupled angular momentum, the matrix 
elements are given by 



In the calculations reported here, we separate neutron 
from proton orbitals, rather than putting them in a single 
chain. As we will soon see, this leads to a three-block 
strategy for implementing the JDMRG truncation. 

With that as background, we now describe how we 
implement the two stages of the analysis, the warmup 
stage and the sweeping stage, for the 48 Cr test problem 
we have considered. 
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type of particle (neutrons, for specificity), building from a 
block B to an enlarged block B' . The medium in which 
we carry out the truncation of that enlarged block in- 
cludes two components. One component, denoted M\ 
in the figure, involves the rest of the neutron orbits, for 
which the corresponding optimal structure was defined 
in the previous sweep (or in the warmup phase). The 
second component, denoted M 2 in the figure, involves all 
of the orbits of the protons, for which likewise the opti- 
mal structure was obtained in the previous stage. Put 
another way, the states in the superblock have the form 



FIG. 3: Schematic illustration of the three-block sweeping 
strategy used in the calculations reported herein. The en- 
larged neutron block B' is immersed in a medium consisting 
of the remaining neutron levels (block Mi) and all the proton 
levels (block M 2 ). 



1. Warmup 

As noted earlier, the initial warmup stage of the 
DMRG procedure involves choosing a first approxima- 
tion to the optimal truncated structure of each group of 
single-particle levels. Because of our separation of neu- 
trons from protons, the groups of relevance are groups 
of neutron orbits and separately groups of proton or- 
bits. Labelling the four orbits as 1 — > 4 (note that we 
have not yet defined the order of the four active orbits), 
we need information on the following blocks: B(l,2), 
5(1,3), -B(l,4). While we need them for both neutrons 
and protons, the symmetry of the problem allows us to 
choose them the same for both. 

Our procedure is to gradually build the identical or- 
bit block, by systematically adding orbit i + 1 to the 
block B(l,i). At each stage of the warmup growth, we 
construct a superblock by coupling the enlarged neutron 
block to the corresponding enlarged proton block (i.e., 
the block with the same set of orbits). We construct the 
superblock for N=Z, from 0^4, and for J v = + . Tar- 
geting the ground state for each of these N=Z systems, we 
obtain the corresponding reduced density matrices and 
implement our truncation. 

The truncation is carried out to those m states with 
the largest density matrix eigenvalues. We do not worry 
about making sure that at least one state is kept for each 
N(or Z) and J, as our experience in these calculations 
indicates it is not necessary. 



2. Sweeping 

As noted earlier, our procedure of separating neutron 
from proton orbitals, rather than lining them up in a 
single chain, gives rise to a three-block strategy in the 
sweeping process. We now expand on what is meant by 
this, by reference to figure 3. 

Let us assume that we are growing the system for one 



i\(n v a v J v ) (4 - n v d v J' v )\ J * (4 & J V )} L ° 

L J M=0 



with n v neutrons in the enlarged neutron block B' , 4— n v 
neutrons in M\ and 4 protons in M^. 

By diagonalizing the hamiltonian in that basis, we ar- 
rive at a ground state wave function 

\L = G.S.) = X^j^j/p^j^ 

a„J„a' v JlP 1 ,J„K„ 

L=0 



X 



\[(n v a v J v ) (4 - n v ol v J' v )\ Ji ' (4 n J n )} 



M=0 



(11) 



from which the reduced density matrix associated with 
B' can be readily evaluated and then diagonalized for 
each n v and J v . Truncation to the m states with 
the largest density matrix eigenvalues gives the optimal 
structure of the truncated block B' , which is stored for 
use in the next sweep stage. 



3. The whole process 

We now have in hand all that is needed to fully de- 
scribe our calculational algorithm. It involves the follow- 
ing steps. 

1. Order the single-particle states. Indeed, in the re- 
sults we will report, we simply use the ordering 
P\/2, P3/2, fe/2, ./V/2 in terms of increasing angular 
momentum, with no effort at optimization. 

2. For each m, we carry out the warmup as prescribed 
in Sec. IV A storing the results (i.e., the reduced 
matrix elements) for each group of orbits. 

3. We then begin by sweeping downwards through the 
orbits, growing from orbit 4 to a block S(4, 3). The 
medium for implementing this truncation involves 
the block -8(1, 2) for the same type of particle and 
the block B(l, 4) for the other, both of which were 
obtained either in the warmup or in the previous 
sweep. 

4. We iterate the above step until all sets of orbits 
have been treated. 
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5. We then reverse and begin sweeping upwards, grow- 
ing the block and implementing a truncation based 
on its coupling to the other two relevant blocks. 

6. We view a sweep down and the following sweep up 
as a single sweeping stage. It involves six growing 
steps, each one giving a ground state eigenvalue. 
We denote the lowest of these six eigenvalues by 
Emini where n denotes the sweeping stage number. 

7. As we systematically sweep up and down, we com- 
pare E^ in with E^ n , *- e > we & t the change 
in the minimum energy from one sweeping stage to 
the next. When this change is sufficiently small, 
typically of order 10 -4 MeV, we stop the calcula- 
tion. No more than 7 sweeping stages were required 
to achieve this level of convergence in any of the 
calculations we carried out. 

8. We then increment to and repeat the set of steps 
2^7. 



V. RESULTS 

Our results for the ground state are presented in table 
1. The exact calculation produces a ground state energy 
of -32.95 MeV. The DMRG calculation with m = 50 
produces a result of —32.80 MeV, about 150 keV from 
the exact result. In this case, the largest matrix that had 
to be diagonalized had a dimension of 3,657. The con- 
clusion is that the three-block JDMRG algorithm gives 
a good reproduction of the exact results for the ground 
state with a fairly small number of states retained in each 
block, but by calculating matrices with a substantial frac- 
tion (roughly 8-9%) of the full basis. 

TABLE I: Results for the ground-state energy from the 
JDMRG calculations described in the text. Max Dim refers 
to the maximum dimension of the hamiltonian matrix that 
required diagonalization. 



m 


Eos (MeV) 


Max Dim 


20 


-32.28 


1,591 


25 


-32.47 


1,893 


30 


-32.57 


2,327 


35 


-32.66 


2,685 


40 


-32.72 


3,109 


45 


-32.76 


3,403 


50 


-32.80 


3,657 


Exact 


-32.95 


41,355 



For any value of to, we can calculate excitation energies 
as well, by working in spaces that are based on the opti- 
mized blocks. As a reminder, during the sweep procedure 
we diagonalize the superblock hamiltonian in a basis in 
which all states of the added level are included. Following 
the density-matrix-based truncation, the enlarged block 
has fewer states, and often significantly fewer. It is in this 



reduced space that we can recalculate the ground state, 
by diagonalizing in the truncated + space, and also carry 
out diagonalizations for other angular momenta. The re- 
sults for to = 50 are illustrated in table 2, for the exci- 
tation energies of the lowest 2+, 4+ and 6 + states. The 
result for the lowest 2 + state is in quite good agreement 
with the exact result, despite the diagonalization of a 
matrix that is less than 5% of the full matrix. For the 4 + 
and 6 + states, the results remain fairly good, although 
progressively worse than for the lower J values. These 
results can of course be improved by further increasing 
to. Overall, we conclude that the method seems capa- 
ble of producing reasonable results for low-lying excited 
states as well, even when we target the ground state only. 

TABLE II: Results for the excitation energies of low-lying 
excited states from the DMRG calculations described in the 
text at a value of m=50. All results are in MeV . 



r 


JDMRG 


Exact 


2+ 


0.84 


0.81 


4+ 


2.02 


1.82 


6+ 


3.95 


3.40 



VI. SUMMARY AND OUTLOOK 

In this paper, we have described an approach for solv- 
ing the nuclear shell-model problem in cases in which 
it cannot be treated by exact diagonalization. The ap- 
proach makes use of the Density Matrix Renormalization 
Group to systematically truncate the space in which di- 
agonalizations are carried out on the basis of well-defined 
dynamical considerations. The approach we develop pre- 
serves rotational symmetry throughout, to avoid losses 
of information associated with breaking of that critical 
nuclear symmetry during the truncation process. 

We test these ideas in the context of a large-scale, but 
exactly diagonalizable, shell-model problem. Namely we 
apply it to the well-deformed nucleus 48 CY with a realis- 
tic KB3 hamiltonian. For these test calculations, we used 
a three-block DMRG strategy, in which neutron and pro- 
ton levels were not mixed in blocks. This leads to certain 
simplifications in the formalism, but at the cost of having 
to treat larger superblock matrices. 

The calculations are able to produce a fairly accu- 
rate reproduction of the exact results for 48 Cr, both for 
the ground state and the first excited 2 + state, both 
members of the ground state rotational band. It does 
this by retaining a fairly small number of states in each 
block. However, because of the three-block approach that 
was used, the dimensions of the superblock matrices we 
needed to treat were still a sizable fraction (5-10%) of 
the full shell-model space. Results for the lowest states of 
higher spins get progressively worse as the spin increases. 

Our work in the future will focus on two fronts. On the 
one hand, we will continue to work with the three-block 
approach, addressing such issues as (i) the possibility of 



7 



improving convergence by altering the order in which lev- 
els are included, building on the work of Legeza and col- 
laborators [Tl|, and (ii) the possibility of speeding up the 
calculations by replacing m as the convergence variable 
by the sum of neglected density matrix eigenvalues [l2| . 
At the same time, we are now in the process of developing 
a two-block JDMRG code 01, in which blocks contain- 
ing both neutron and proton orbits are built. This should 
enable us to dramatically increase the number of states 
retained in a block for a given size of the full superblock 
matrix. This will be critical as we consider the applica- 
tion of these methods to heavier nuclei, the ultimate goal 
of the project. 
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